The Compressible Ising Spin Glass: Simulation Results 
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This paper reports numerical studies of a compressible version of the Ising spin glass in two 
dimensions. Compressibility is introduced by adding a term that couples the spin-spin interactions 
and local lattice deformations to the standard Edwards- Anderson model. The relative strength of 
this coupling is controlled by a single dimensionless parameter, fj,. The timescale associated with the 
dynamics of the system grows exponentially as [i is increased, and the energy of the compressible 
system is shifted downward by an amount proportional to fj, times the square of the uncoupled 
energy. This result leads to the formulation of a simplified model that depends solely on spin 
variables; analysis and numerical simulations of the simplified model predict a critical value of the 
coupling strength above which the spin-glass transition cannot exist at any temperature. 

PACS numbers: 75.10.Nr,75.40.Mg,05.50.+q 



I. INTRODUCTION 

Much theoretical study has been made of the nature 
of the spin-glass transition. It is now generally accepted 
that the three-dimesional spin glass undergoes a second- 
order phase transition at finite temperaturejiiSiiiiiSiSi! 
and the bulk of the evidence in two dimensions is consis- 
tent with a zero-temperature phase transition^&^iiSiiiii^ 
although recent work suggests that the lower critical di- 
mension for some spin-glass models is greater than two^» 
The continued controversy hints at the delicate and sub- 
tle nature of the spin-glass transition and suggests that 
modifications to the underlying model, even small ones, 
could have dramatic effects on the system. 

Compressibility has already been shown to have a 
strong effect on a variety of spin systems. The inclu- 
sion of compressibility in the Ising ferromagnet modi- 
fies the standard second-order transition to a first-order 
transition that occurs at the Curie temperatureii^^ The 
(fully frustrated) 2-D triangular Ising anti-ferromagnet 
does not undergo a phase transition; however, when 
compressibility is added to the model, the character- 
istic frustration is relieved, and the system develops 
a first-order transition to a "striped" phase at low 



temperatures 



16.17.18 



Other frustrated spin systems are 



known to have their frustration relieved by the presence 
of magnetoelastic couplingSfiS^ and polaron effects alter 
the nature of magnetic transitions in frustrated physical 
systems such as manganitesiSii 2 ^ 2 ^ 

With the possibility of relieving frustration, the addi- 
tion of compressibility to spin-glass models could dramat- 
ically alter the nature of the spin-glass phase and/or the 
transition thereto. Furthermore, the fact that all physi- 
cal systems must possess some (albeit small) spin-lattice 
coupling provides a physical motivation for such studies. 

A previous paper introduced a particular model for 
the compressible spin glass with a linear coupling be- 
tween the spin-spin interactions and the distances be- 
tween neighboring particles^* The work described there 
involved simulations of the compressible spin glass per- 
formed on two-dimensional systems in which the volume 



was held fixed. Results of the direct simulations sug- 
gest a simplified model, qualitatively equivalent to the 
first, that depends only upon spin degrees of freedom. 
The presence of compressibility alters the preferred spin 
configurations of the system, so that the transition to 
a low-temperature spin-glass phase is impossible above a 
critical value of the coupling. The current paper expands 
on that previous work as well as provides details of the 
analysis. Presented here are results showing the expo- 
nential slowing down of the time to reach equilibrium 
as the coupling increases, additional quantitative moti- 
vation for the simplified model, and a functional form 
for the entropy of the spin glass, from which thermody- 
namic quantities are predicted. Finally, a phase diagram 
illustrates an approximate boundary separating critical 
behavior from the region where the spin-glass transition 
cannot exist. 

The structure of this paper is as follows: Section [H] 
describes the Hamiltonian of the compressible spin glass 
and defines the important tuning parameters. The details 
of the computer simulations are discussed in Scc. lIIII and 
the results of the simulations are presented in Sec. IIVI 
A simplified model is introduced in Sec. [V] along with 
results of numeric simulations and analytic investigations 
performed on the simplified model. Section IVII contains 
the primary conclusions and some additional points of 
discussion. 



II. THE MODEL 



The Hamiltonian for the compressible Ising spin glass 
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ft = ~ JijSiSj + a ^2 JijSiSj (rij - r ) + [/lattice ■ 

(i.j) (id) 

(1) 

The first term is the standard Edwards- Anderson spin- 
glass Hamiltonian^ with the sum performed over pairs 
of nearest neighbors. The spins Si are dynamic variables 
which may take the values +1 or — 1. The interactions Jy 
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are chosen randomly from {± J} with equal probability 
and are then held fixed; this collection of interactions 
represents a single realization of the quenched disorder 
central to the nature of the spin glass. 

The coupling between the spin interactions and the 
lattice distortions is contained within the second term 
of Eq. lJT]l. where the coupling is considered to linear 
order with proportionality constant a. This constant 
multiplies the change in bond length: r»j represents the 
Euclidean distance between particles i and j, and r is 
the natural spacing of nearest neighbors on the lattice. 
This term allows the system to lower the total energy 
by displacing the particles from their regular lattice po- 
sitions. Spins with satisfied interactions (i.e., those with 
JijSiSj = +1) will tend to move closer together in order 
to strengthen the effect; similarly, unsatisfied bonds will 
tend to lengthen as the particles move farther apart to di- 
minish the negative effect of their interaction on the total 
energy. The inability of all of the bonds in the system to 
distort simultaneously in the ideal fashion is the mech- 
anism by which the degeneracy of configurations with 
equal spin-spin energy is broken. 

The final term, [/lattice, stabilizes the lattice by pro- 
viding a restoring force to counteract the displacements 
generated by the spin-lattice interactions. The stabiliza- 
tion is obtained by connecting harmonic springs between 
nearest neighbors and between next-nearest neighbors 
(along the diagonals of the square lattice). Each spring 
has as its unstretched length the natural spacing of the 
vertices, so that Lattice is zero in the absence of spin- 
lattice coupling when the particles are not displaced. 

Two important parameters can be formed. The first 
arises due to force balance between the last two terms in 
Eq. (TJ): 



This has dimensions of length and represents the scale 
of the typical displacements of the particles from their 
uncoupled locations on the square lattice. The second 
parameter is 



which is dimcnsionless; it represents the strength of the 
spin-lattice coupling, relative to the spin-spin interaction. 
When fi — 0, there is no spin-lattice coupling, lattice 
distortions are not energetically favorable, and the model 
reduces to the standard Edwards-Anderson spin glass. 
The interaction strength J serves merely to set the energy 
scale for the model. In this work, J is set to unity, while a 
and k are chosen so that 5 and fi take the desired values. 

III. SIMULATION DETAILS 

All simulations were run on square lattices of linear 
dimension L with periodic boundary conditions in two 



dimensions; the size of the system was held constant in 
each direction, fixing the total volume. The control pa- 
rameters were adjusted so that S was set at ten percent 
of the natural lattice spacing while \i was varied over 
the range < fi < 5. Since the information regarding 
relative energy scales is contained within /i, the specific 
value of S does not affect the qualitative nature of the 
results, so long as the displacements are small enough to 
maintain the topology of the lattice. As discussed below, 
however, small nonlinearitics do depend on the extent of 
the lattice distortions. 

For all values of the control parameters, 100 different 
bond configurations (i.e., realizations of the quenched 
disorder) were simulated. Calculated quantities were 
then averaged over the various runs. 

Two different methods of simulation were used to study 
the compressible spin glass. For the first method, suitable 
for studying the dynamics of the model, states were gen- 
erated via single-spinflip Monte Carlo steps, with transi- 
tion probabilities dependent upon the difference in energy 
between the two spin states. These energies employed the 
full Hamiltonian of Eq. including the components 
that depend on the particle positions. For purposes of 
determining transition probabilities, the spins were con- 
sidered to flip in place, i.e., without any particle motion. 
The lattice was then relaxed for the new spin configu- 
ration. The full simulation algorithm is as follows: The 
system is started in a random spin configuration with 
the particles located at the positions which minimize the 
total energy. From a given spin configuration, a particle 
is chosen at random. This particle is given a chance to 
flip, in place, from the state with energy E\ to the state 
with energy E 2 - If E 2 < E\ the spin is flipped; other- 
wise the spin flips with probability exp[— (E 2 — E\)/T]. 
After L 2 randomly chosen particles have been considered 
(i.e., one Monte Carlo step), the lattice is relaxed to the 
minimum of the potential energy for the new spin config- 
uration using conjugate-gradient minimization. System 
properties are recorded for analysis, and this process is 
then repeated. 

In order to ensure proper equilibration, I follow the 
algorithm prescribed by Bhatt and Youngi^ The spin- 
glass susceptibility 

Xsg = Yp. i^iSj) , 

where (•) indicates a thermal (time) average, is calculated 
by two different methods, each as a function of equilibra- 
tion time iequii- One method uses the overlap between 
states of the same system at two different times during 
the run, while the other uses the overlap between states 
of two randomly initialized, independently run replicas 
of the same bond realization. These two computation 
methods produce the same value of x sg as £ oqu ii — » oo, 
but the "two-times" method approaches the asymptotic 
value from above, while the "two-replicas" method ap- 
proaches from below. When the two values are within 
statistical error of one another, the system is equilibrated. 
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The simulations are typically run for several multiples of 
the equilibration time in order to acquire data from un- 
corrected portions of the time evolution. 

Another method of simulation, suitable for studying 
static properties such as the energy, involves substitut- 
ing a collection of pre-generated spin states into the 
compressible spin-glass Hamiltonian and relaxing each 
to the minimum of its total energy with respect to the 
particle positions. Typically, the spin states are gener- 
ated by single-spinflip Monte Carlo simulations using the 
standard (incompressible) spin-glass Hamiltonian, which 
takes much less time than simulations of the full Hamilto- 
nian as described above. In this manner, "typical" states 
may be analyzed to observe the effect of the compressible 
terms on quantities of interest; however, these states will 
not occur with frequency given by the correct Boltzmann 
weight, so care must be taken not to draw conclusions 
that would rely on such an assumption. For the small- 
est system sizes (L — 3, 4, and 5), it was possible to 

T 2 

enumerate all 2 possible spin states for a given bond 
configuration. 

For all methods of generating spin states, the lattice 
was relaxed to its minimum using the conjugate-gradient 
minimization technique. 26 Since the distortions of the 
lattice are kept small by the value of 5, the potential- 
energy landscape is close to quadratic, and the minimum 
can typically be located to reasonable numerical toler- 
ance within a few conjugate gradient steps. Neverthe- 
less, because of the computation involved in calculating 
the lattice energy, this portion of the simulation takes 
approximately two orders of magnitude more time than 
the Monte Carlo spinflips. 



IV. RESULTS 

Simulations of the two-dimensional, constant-volume 
compressible Ising spin glass were performed for system 
sizes ranging from L = 3 to 40 using the techniques de- 
scribed above. Data from these simulations are presented 
and anlayzed below. 



A. Dynamics 

The time required for the system to reach thermal equi- 
librium is an easily accessible measure of the timescale 
for the system dynamics. For each value of /z, a dif- 
ferent equilibration time is required, and Fig. ^ shows 
the dependence of the equilibration time, t cqu ii, on [i for 
the L — 10 systems at a relatively high temperature, 
T = 2.0. As the fit line on the semilog plot demon- 
strates, the growth of the equilibration time, in Monte 
Carlo steps (MCS), is exponential in /i; the slope of the 
exponential fit is 1.8 MCS -1 . The rapid growth of the 
equilibration time as the coupling is increased can be 
viewed as a growth of energy barriers between states that 
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FIG. 1: The time to reach equilibrium, i e quii, grows expo- 
nentially as increases. The data here are from 100 L = 10 
systems at T = 2.0, and the slope of the exponential fit line is 
1.8 MCS -1 . The dramatic increase of simulation time makes 
straightforward simulation of the dynamics difficult for large 
values of the coupling. 

were previously similar in energy. The movement of par- 
ticles "locks in" the current spin configuration, increasing 
the timescale for single-spinflip transitions. 

The growth of the equilibration time as a function of 
the coupling is in addition to the usual dramatic growth 
of dynamic timescales as the temperature is lowered (see, 
for example, Fig. 3 of Ref. 0). Since the number of simu- 
lation steps required increases exponentially with /i and 
the computation time per step increases in a manner pro- 
portional to the number of spins, direct simulations of the 
system dynamics at temperatures approaching the tran- 
sition become prohibitive for large values of the coupling. 

B. Energy Analysis 

In analyzing the results of the simulations, the vari- 
ous components of the total energy may be computed 
independently for a given spin configuration. Of partic- 
ular interest is the first term in Eq. . This component 
represents the contribution due solely to spin-spin inter- 
actions and is denoted Eq. It is equivalent to the energy 
of that spin configuration on an undistorted lattice in the 
absence of any spin-lattice coup ling. 

As shown in Fig. 1 of Ref. [2J, the effect of the coupling 
is to shift the states of the system downward in energy. 
When [i = 0, the energy levels are 5-functions separated 
by constant gaps of 4J, the smallest energy difference 
between states of the incompressible ±J model. As fi is 
increased from zero, each energy level (identified by Eq) 
shifts downward in energy and broadens into a Gaussian- 
shaped band. 

For all of the states with a given value of Eq , the dis- 
tribution of energies is characterized by two values: the 
average shift in energy, AE(Eo,fi) = (E(Eo,/j,)) — Eq, 




FIG. 2: The effect of the coupling on the energy for a single 
fully enumerated L — 4 system, (a) The average shift in 
energy, AE, is plotted as a function of the coupling ji. For 
each band of states, from Eq = — 24 (the ground-state energy 
for this particular system) to Eq = 0, the energy shifts by an 
amount proportional to the coupling, (b) The width of each 
band, a, is also proportional to /x. 



and the width <j(Eq) as given by the standard deviation 
of the distribution. Both AE and a are linearly propor- 
tional to /i, as shown in Fig. The data in that figure 
were obtained from a single L — 4 system using com- 
plete enumeration of all spin configurations; each line 
represents the data for a value of Eq ranging from the 
ground-state energy for this specific system, E = —24, 
to Eq = 0, where there are equal numbers of satisfied and 
unsatisfied bonds. 

The proportional dependences of both the energy shift 
and the width on /i are due to the fact that each spin state 
individually shifts by an amount exactly proportional to 
the coupling. When minimizing the potential energy of 
the lattice for a given spin configuration, the positions of 
the particles are determined by the value of 5; the value 
of fi then multiplies the result to determine the total 
energy in the distortions. Due to this fact, it is possible 
to characterize changes to the energy of the system at 
any convenient value of [i and then scale the obtained 
quantites by the coupling. 

The lines shown in Fig. [3 have different slopes, indicat- 
ing that the various bands shift and broaden at different 
rates as n increases. The states with higher Eq move 
downward in energy more rapidly than lower-i?o states. 
Data for the shift in average energy from the uncoupled 
value, scaled by [A, are plotted as a function of original 
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FIG. 3: Dependence of the energy shift on the spin-spin en- 
ergy Eq. These data are averaged over 100 L = 10 systems 
run at a variety of temperatures. The parabolic shape of this 
curve results from the fact that configurations with roughly 
equal numbers of short (satisfied) and long (unsatisfied) bonds 
can distort more effectively than those with many bonds of 
the same length. 

energy level, Eq, in Fig. 100 systems with L — 10 and 
H = 0.1 were run at a sequence of temperatures and av- 
eraged to produce this plot. In practice, the states with 
positive Eq are difficult to populate at finite temperature 
due to the exponential suppression of the Boltzmann fac- 
tor. 

The parabolic form of this curve can be explained 
by the observation that with the volume held con- 
stant, configurations with predominantly short (or long) 
bonds cannot distort as effectively as configurations with 
roughly equal numbers of short and long bonds. For even- 
valued system size L, this curve should be symmetric 
about Eq = since there is a relationship between spin 
states with alternate spins flipped: long bonds become 
short bonds and vice versa, resulting in a state with Eq 
of equal magnitude but opposite sign that has an identi- 
cal energy shift. 

The lack of exact symmetry about Eq = is due to 
small nonlinearities resulting from non-zero 5. Figure 0] 
shows data for the typical value of S = 0.1 along with a 
sample of data in which 6 was set to 0.01. The results are 
qualitatively similar, though the smaller-distortion curve 
is more symmetric. 



C. Size Dependence 

To study the size dependence of the energy, data for 
AE and a was collected for system sizes from L = 3 
to 40. For all system sizes, the shift in the average en- 
ergy displays the parabolic shape shown in Figs.|3|and21 
and the similarity in form suggests that the curves may 
be made to collapse. Fig. 02a) contains the results for 
AE I [i for the full range of system sizes simulated. 100 
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FIG. 4: The effect of changing the distortion parameter 8, as 
defined in Eq. @, is shown for the fully enumerated systems 
with L — 4 and /j, — 0.1. The parabolic form of the data is un- 
changed; however, the smaller value of the typical distortion 
size results in a curve that is slightly more symmetric about 
Eo ~ 0. 100 systems were averaged to produce these data. 



systems were averaged at each system size; in the figure, 
data points are only displayed for values of Eq where at 
least 20 systems were represented. The L — 3, 4, and 5 
systems were fully enumerated, while the larger systems 
were run at a series of temperatures to obtain data over a 
range of values of Eq. As the inset in that figure demon- 
strates, when both axes are scaled by L 2 , the data for the 
various system sizes approach a constant curve as L in- 
creases. While there are finite-size effects in the smallest 
systems, the data for L > 10 collapse quite well. 

The quadratic form of the scaled data for AE is ex- 
pressed as 



AE 

JIl 2 



= Aae x 



Eo 
L 2 



- B 



AE 



a 



AE ■ 



(4) 



The locations of the minima for each system size were av- 
eraged to determine the global horizontal offset: Bae = 
0.063 ± 0.001. With Bae determined, the L = 40 
data were then fit to the parabolic form above, with 
Aae = 0.1166 ± 0.0007 and Cae = -0.5004 ± 0.0008. 
The main panel of Fig. Ufa) shows the data and fit plot- 
ted in a manner that makes the collapse to the form of 
Eq. |0J apparent. 

Data for the width of each band also demonstrate a 
quadratic function of Eq , as shown in Fig. \5j[b) . As with 
the energy shift, the a data for the various system sizes 
can be scaled to lie on a common curve; however, while 
the Eq axis is again scaled by the system size L 2 , the 
width axis is only scaled by the linear size of the system, 
L. 

The scaled data for the width are described by the form 



c n 



(5) 
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FIG. 5: System-size scaling, (a) The slope of the average 
energy shift, as a function of Eo, with both axes scaled by L 2 . 
As L increases, the data for different system sizes collapse 
onto a common curve that is quadratic in Eo- The main 
panel shows the data plotted and fit according to the form of 
Eq. @, while the inset shows the scaled data directly, (b) The 
data for the spread in energy as a function of Eq can also be 
made to approach a common parabolic curve; however, while 
the Bo-axis is again scaled by L 2 , the width axis is scaled by 
the linear size only. 



As with AE, data from all sizes were used to obtain B a = 
0.039 ± 0.004. The L = 40 data were then fit to Eq. ®, 
resulting in A a = -0.085 ±0.002 and C a = 0.308 ±0.003. 
Figure E^b) shows the data and fit line. 

The scaling behavior of a implies an interesting side 
effect of the introduction of compressibility. Since Eq 
is proportional to L 2 , the total number of spins, the 
right-hand side of Eq. JSJ is independent of L, and thus 
a ~ jiL. Neighboring energy bands will overlap to a large 
degree when the width of the bands is comparable to the 
spacing between them, i.e., when fiL ^4. As L — > oo, an 
infinitesimal value of the coupling will satisfy this condi- 
tion, rendering the previously discrete energy spectrum 
continuous. 

While the forms of AE and a are similar, it is not 
immediately apparent that the two quantities are directly 
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FIG. 6: The ratio of the width to the magnitude of the energy 
shift as a function of Eo for different system sizes; the axes 
have been scaled by L and L 2 , respectively, to demonstrate 
an approach to constant behavior as the system size increases. 
Plotted this way, the scaled ratio is less than a constant of 
order unity. Thus, a becomes negligible compared to AE for 
large L. 

related. In fact, over a large range of Eq, the spread in 
energy is proportional to the energy shift, as shown in 
Fig. El where the data displayed in Fig. are plotted as 
a ratio of a x L to the absolute value of AE versus E /L 2 . 
Again, the data from different system sizes were made to 
collapse by appropriate scaling of the axes. 

It is apparent from Fig. that the magnitude of the 
scaled ratio is less than a constant, p, with p rj 0.6 for 
large values of L. The relationship between a and AE 
can be expressed as 

o- < £ |AB| . (6) 

Thus, as the system size increases, the width of a band 
of states becomes negligible compared to the magnitude 
of the shift in energy from its uncoupled value. 

V. SIMPLIFIED MODEL 

The form of AE, as demonstrated in Figs. E3 El and(5{a) 

and expressed by Eq. Q, is a quadratic function of Eq. 
In addition, the spread in the energy becomes negligible 
compared to the energy shift for large system sizes, as 
Eq. © demonstrates. These observations motivate^ an 
approximate Hamiltonian for the compressible spin glass: 

Wapprox = — J]] JijSiSj + — I JijSiSj J . (7) 

The sum that appears in both terms is performed as de- 
scribed for the original Hamiltonian of Eq. Q , and nu- 
merical factors — such as j4ae from Eq. (0} — have been 




FIG. 7: A selection of data from various system sizes shows 
the relationship between the two position-dependent compo- 
nents of the total energy given by Eq. here, the second 
term (representing the energy due to the coupling, -E C ou P iing) 
is plotted versus the third term (the lattice energy, -Biattice) 
for individual spin configurations. Both axes are scaled by L 
to bring the data from different system sizes into a common 
range. As the coupling energy is proportional to the lattice 
energy, these two terms may be combined into a single term 
that describes the energy shift due to particle motions in the 
presence of the coupling. 



absorbed into the coupling constant so that v w 0.12/^. 
Again, the first term represents the energy due to spin- 
spin interactions, denoted Eq. The second term, which 
contains the coupling between the spins and the lattice, 
can be viewed as a combination of the final two terms 
of Eq. (JJJ with a typical distortion of the order of 5 as 
defined in Eq. ©. 

That the two position-dependent components of the 
total energy may be combined in this way is shown ex- 
plicitly in Fig. where the coupling energy is plotted 
against the lattice energy for data obtained in the previ- 
ously described simulations. Both axes are scaled by L 2 
to bring the points from different system sizes into a com- 
mon range, and the solid line is a linear fit to all of the 
data. The coupling energy is proportional to the lattice 
energy, and since the two terms are related in a straight- 
forward manner, their net effected may be represented 
by a single term in the simplified Hamiltonian. 

I note some features of the approximate model. First, 
the system size must be included explicitly in order to 
preserve the observed scaling behavior. Second, rather 
than being constructed from a combination of parame- 
ters, the coupling constant v is directly present and con- 
trols the strength of the compressibility term. Finally, 
and most importantly, this Hamiltonian contains only 
spin degrees of freedom; the positional variables are ab- 
sent, and the degrees of freedom associated with them 
have been absorbed into the second term of Eq. J7J. 
Thus, the simplified model can be viewed as a mean-field 
version of the original Hamiltonian, where the energy due 
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FIG. 8: Results from Monte Carlo simulations of the approx- 
imate model of Eq. J7J on 100 systems with L — 10. The 
average value of Eo is plotted as a function of T for various 
values of v. Below a critical value, v* ~ 0.35, the energy 
approaches that of the ground state (dashed line) as T — » 0; 
above u* , the energy limit is predicted by Eq. IjlOfl . The 
solid lines are not fits; rather, they represent predictions ob- 
tained by minimizing the free energy using the entropy form 
of Eq. 1141 . This figure is reprinted from Ref. I24I 

to local distortions has been replaced by an energy con- 
tribution that is typical for states with equivalent spin 
energy. As a practical matter, this feature also means 
that simulations and analytical work may be performed 
on the model using techniques identical to those used in 
standard (incompressible) spin-glass studies. 



A. Simulation Results 

To simulate the model described in the simplified 
Hamiltonian, singlc-spinflip Monte Carlo simulations 
were performed at various values of the coupling as 
parametrized by v. As in standard Monte Carlo spin- 
glass simulations, for each bond realization a sequence 
of spin states was generated at a fixed temperature T, 
with transition probabilities between states based on the 
difference in energy as calculated from Eq. (7J. 

To monitor which states of the compressible spin glass 
are favored at a given value of T and v, the thermal- 
and disorder-averaged values of E , denoted (Eo), are 
calculated. Results for 100 systems with L — 10 are 
shown as data points in Fig. [H] Solid lines in that figure 
represent predictions based on the free-energy analysis 
described below in Sec. IV CI 

For each value of u, the average value of Eq decreases 
as the temperature is lowered (at infinite temperature, 
(Eq) = since the Eq = states have the highest en- 
tropy). For small values of v, as T — > 0, the data ap- 
proach the ground-state energy of the uncoupled system, 
indicated by the dashed horizontal line in the figure. 
For each curve as a function of T, the inflection point 



represents the temperature at which the presence of the 
ground states of the uncoupled system becomes impor- 
tant. Below this temperature, those ground states begin 
to be populated, halting the decrease in (Eo). As the 
coupling is increased, the inflection point moves down 
in temperature, eventually disappearing when v w 0.35. 
For larger values of the coupling, there is no inflection 
point, and the data for (Eq) no longer approach the 
ground-state energy as T — > but rather terminate at 
some higher value that increases with increasing v. 

The heat capacity can be calculated as the fluctuations 
in the energy about its average value, divided by the 
square of the temperature. For the simplified model of 
the compressible spin glass, the specific heat shows no 
interesting features, going smoothly to zero as T — > 0. 
However, a similar quantity, using Eq instead of the total 
energy, can be calculated: 



Data for C from the simulations of the simplified model 
are shown as points in Fig.^a). Solid lines in that fig- 
ure are predictions based on the free-energy analysis de- 
scribed below. 

At small values of v, the data and corresponding pre- 
diction for C are peaked. For v = 0, i.e., the standard 
Ising spin glass, this peak is interpreted as signaling the 
onset of critical behavior that precedes the spin-glass 
transition as the temperature continues to be lowered^ 
As the coupling approaches a critical value, v*, the tem- 
perature at which the peak in C is located moves toward 
zero, and the height of the peak, C max , diverges. Fig- 
ure Efb) shows that this divergence is a power law: 

C max = A(v*-is)- p . 

The points in that figure are the locations of the maxima, 
as obtained from parabolic fits to data near the peak of 
each curve as a function of T, while the solid line is a 
power-law fit to the peaks of the predicted curves with 
v* = 0.365 and a power-law exponent of 1.3. The reasons 
for the divergence are discussed below. 



B. Analytic Results 

For the approximate Hamiltonian of Eq. Q, there is 
a one-to-one correspondence between the spin energy Eq 
(calculated as before) and the total energy. To under- 
stand the results of the simulation as presented in Fig. 03 
it is useful to analyze the expected value of the energy 
for various values of the control parameters. In simplified 
notation, Eq. can be written^ 

E = E + ^Eo 2 . (9) 
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FIG. 9: (a) The heat-capacity-like quantity C — defined in 
Eq. (|HJ — as a function of T for different values of v, calcu- 
lated from simulations of the simplified model on 100 aver- 
aged systems with L — 10. The peak in C shifts downward 
in T as v increases. The solid lines are predictions based on 
the free-energy analysis described in Sec. IV Cl (b) The max- 
imum value of C demonstrates a power-law divergence as v 
approaches a critical value, v* . The points are obtained from 
polynomial interpolation of data near the peak of each curve 
in panel (a). The solid line is a fit to the peaks of the C 
prediction curves as a function of v\ from the fit, v* = 0.365, 
and the power-law exponent is 1.3. 



Taking the derivative with respect to Eq allows one to 
calculate the value of the spin energy, E^mm, that mini- 
mizes the total energy as function of v. 



O.min 



—L 2 



(10) 



This represents the expected value of the spin energy as 
T — > 0; however, the calculated value is not necessar- 
ily realizable, since Eq must be greater than the ground- 
state spin energy E . gn d- For small values of v, the calcu- 
lated value of -Eo im in lies in the non-physical region below 
^o.gnd, and thus the minimum spin energy is, of course, 
equal to the ground-state spin energy. This explains the 
form of the small-j/ data in Fig. [81 where the average 
spin energy decreases with temperature but approaches 
-^o.gnd asymptotically as T — > 0. 

There is a value of v at which the minimum Eq be- 
comes equal to the ground-state energy of the uncoupled 



-i 2 



2E, 



0,gnd 



(11) 



For the two-dimensional spin-glass, the ground-state en- 
ergy per spin is —1.4^ and thus v* = 0.36. This value is 
the same as the value of v at which the simulation data, 
as shown in Figs. |H1 and display a change in behavior. 

At v = v* , the nature of the energy spectrum is dra- 
matically altered: As the coupling increases from zero, 
higher- Eq states shift downward in energy more rapidly 
than lower-En states, and thus the difference in total 
energy between neighboring Eq levels becomes smaller. 
The value v* represents the coupling at which the ground 
state and first excited state of the uncoupled system have 
the same total energy. Above this value, states with 
Eq = -Eo.gnd no longer have the lowest total energy, and 
as v increases still more, increasingly higher En-levels are 
associated with the ground states of the compressible sys- 
tem. The data in Fig. [HI display this feature, as the zero- 
temperature value of (Eq) increases with the coupling 
for v > v* . The divergence in C, as shown in Fig. EI is 
also a consequence of this change in the energy spectrum. 
As v — > v* , the difference in total energy between levels 
near (Eq) goes to zero, and thus the fluctuations in Eq 
no longer vanish as T — > 0. 



C. Free Energy Analysis 

In order to predict which states are preferred as a func- 
tion of v and T, the free energy must be minimized, and 
for this a functional form for the entropy is needed. The 
probability of generating a state of given energy is pro- 
portional to the Boltzmann-weighted density of states: 
P(E) oc Q(E) exp(-E/T); from this the entropy S(E) is 
derived as logfi. Since the density of states is a function 
of the uncoupled energy, i.e., O = Cl (E(E )), it may be 
calculated for any value of the coupling. For a given v 
and T, simulations will produce a limited range of ener- 
gies that will be populated with statistical significance; 
data is therefore acquired at different couplings and tem- 
peratures to produce overlapping regions of data that 
may be combined. Since for each run the proportionality 
between the generated probabilities and the density of 
states is unknown, it is more convenient to generate the 
derivative of the entropy: 



dEQ-dEQ^ P ^ 



E(Eq)/T] 



(12) 



Data for the derivative of the entropy is shown in 
Fig. ^1 For each system size, 100 individual bond con- 
figurations were run at a variety of temperatures and 
averaged. By plotting dS/dE Q versus Eq per spin, the 
data from different system sizes are made to lie on a sin- 
gle curve. This curve is linear over a large region passing 
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FIG. 10: Data for the derivative of the entropy, averaged over 
100 systems each at sizes from L = 4 to 40. When plotted 
as a function of Eo per spin, the data lie on a common curve 
that is well fit by the functional form given in Eq. 1131 . 



through Eo = 0, and the deviations from linearity are ex- 
ponential. As demonstrated by the solid line in Fig. 1101 
the overall curve is well fit by the functional form 



dS 
~dE~o 



"Cl 



Eo 
L 2 



c 2 sinh c 3 



E 
'L 2 



(13) 



where c\ = 0.5, c-i — 4.5 x 
entropy is thus of the form 



io- 



and C3 = 5.6. The 



S (E ) = So - Si£ 2 - S 2 cosh (S 3 E ) 



(14) 



with S± = 0.25/L 2 , S 2 = 8.1 x 10- 5 £ 2 , and S 3 = 5.6/L 2 . 

With the entropy given by Eq. l|14|l and the energy 
given by Eq. ©, the free energy, F = E(E ) - TS{E ), 
may be minimized with respect to Eq- The spin energy 
of the system in the thermodynamic limit is thus given 
by the solution to the equation 



2i/ 
L 2 



Eo + 2TSiE + TS2S3 sinh (S 3 E Q ) = . 



While this equation cannot be solved analytically, it is 
possible to obtain a numerical solution as a function of 
v and T. Such results are plotted as solid lines in Fig. |H1 
where the values of (Eo) as predicted from the free- 
energy calculation are in excellent agreement with those 
obtained from simulations of the simplified model. The 
predictions tend to diverge from the data at low temper- 
atures for small values of the coupling; this is the regime 
in which the lowest-energy states are heavily populated 
and the functional form for the entropy — which contains 
no low-energy cutoff — ceases to be a good description of 
any actual system. 

The free energy described above was also used to pre- 
dict values for the heat-capacity-like quantity C, as de- 
fined in Eq. (JSJ) . These predictions are shown for various 
values of the coupling as the solid lines in Fig. Efa); the 



0.0 



0.1 



0.2 



0.3 



0.4 



0.5 



0.0 



FIG. 11: Phase diagram in the v — T plane showing the tem- 
perature at which the predicted energy crosses the ground- 
state energy and the temperature at which the peak in C 
occurs. Dashed lines indicate linear fits to these data; the 
two lines terminate at v* . The region beneath the lines rep- 
resents the critical regime, signaling the onset of the spin-glass 
phase. At the same value v* , the order parameter (f> — defined 
in Eq. 115t — increases linearly from zero, indicating the sup- 
pression of critical behavior. 



agreement with data from the simulations is very good 
except at low temperatures for small where the appli- 
cability of this form for free energy is expected to break 
down. Near the peaks in C, however, the predictions 
are still in reasonable agreement with the simulation re- 
sults. The peaks of the predicted curves were thus used 
to generate a prediction for the divergence of C max - The 
power-law fit to these peaks is shown as the solid line 
in Fig. 0Hb), where the power-law exponent is 1.3 and 
v* = 0.365, consistent with the value of 0.36 obtained 
via Eq. lfTT|) . 



D. Phase Diagram 

It is possible to map out various quantities as a func- 
tion of v and T . Consider the temperature at which the 
predicted value of (E ) crosses the ground-state energy 
of the uncoupled system (see Fig.|SJ). This represents the 
breakdown of the prediction due to the lack of a consis- 
tent analytic form for the entropy near the ground state. 
Near this temperature, there is an inflection point in the 
curve of the average Eo as the presence of the ground 
state becomes important and the curve begins to flatten. 
Also of interest is the location of the peak in C with re- 
spect to T, which indicates the onset of critical spin-glass 
behavior. Both of these quantities, calculated from the 
predictions described in Sec. IV CI are plotted in a v — T 
phase diagram in Fig. 1111 with linear fits to the points. 

The lines in the v— T plane mark an approximate phase 
boundary between the normal (paramagnetic) phase and 
the critical regime that signals the onset of spin-glass be- 
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havior. The fit lines for the ground-state crossing and the 
peak in C terminate at v = 0.359 and 0.353, respectively, 
consistent with the value of v* predicted from Eq. 1)11(1. 
Commensurate with the termination of these lines at v* 
is the growth from zero of an order parameter character- 
ized by the difference between the minimum value of Eq 
and the ground-state energy of the uncoupled system: 

4> = £umin — -Eo.gnd ■ (15) 

Using Eqs. i|10|) and (|llfl . it is apparent that 4> grows as 
\/i>* — i.e., linear just above the critical value v* . 
Figure 1771 shows this phenomenon. 

The interpretation of <f> is as a measure of the inac- 
cessibility of low-Eh states, even at low temperatures, 
due to the presence of the coupling to lattice distortions. 
These states, at and near the ground state of the un- 
coupled system, are no longer the lowest-energy states of 
the compressible spin glass, and the competition between 
energy and entropy no longer exists. Non-zero <fi is thus 
correlated with the suppression of critical behavior that 
precedes the spin-glass phase; above v* , the spin-glass 
transition cannot exist. 



VI. CONCLUSIONS 

This paper has analyzed a model^i for a compressible 
Ising spin glass that lends itself to simulations similar 
to those for standard spin glasses, with additional steps 
to determine the positions of the spin particles. While 
exploration of the dynamics of this system has proven 
difficult, it is possible to characterize the effect of the 
coupling to lattice distortions on the energy of the sys- 
tem. Both the shift in energy and the width of each band 
of states display parabolic shapes as a function of the un- 
coupled energy Eq , and system-size scaling demonstrates 
that the width becomes negligible compared to the en- 
ergy shift as L increases. 

The form of the shift in energy due to the presence 
of compressibility motivates a simplified model for the 
compressible spin glassi^i This model, which depends 
only upon spin degrees of freedom, was simulated using 
standard techniques. In addition, analysis of the simpli- 
fied model suggests a critical value of the coupling above 
which the nature of the energy levels changes dramati- 
cally, and the simulation data confirm this. Due to the 



elimination of the critical regime, a spin-glass transition 
cannot exist above this critical value. 

The simplified model adds long-range interactions to 
the nearest-neighbor behavior of the standard Edwards- 
Anderson model for the spin glass. This provides a con- 
venient mechanism for incorporating phonon effects into 
theoretical spin-glass studies, and it is possible that con- 
sideration of these effects may help to shed light on some 
experimental results that have yet to be fully explainedr. 
For example, work by Bitko et aliSl demonstrated the 
existence of a signature for the spin-glass transition at 
high frequencies, hinting that coupling to high-frequency 
modes (such as phonons) may be important. 

I expect the results not to change qualitatively in three 
dimensions. Preliminary studies similar to those de- 
scribed above suggest that, as for the two-dimensional 
case, the shift in energy is proportional to the coupling, 
/j, and to Eq 2 . Furthermore, the energy shift scales as 
the volume of the system, L 3 , while the spread in each 
energy band scales as L. Thus, the assumptions that led 
to the simplified model of Eq. hold even more strongly 
in three dimensions; similar results for the elimination of 
the spin-glass transition above a critical value of the cou- 
pling should then follow, but additional work is required 
to verify this. 

It is important to note that many of these results are 
expected to be quite different if the constraint of con- 
stant volume is removed. The quadratic nature of the 
energy shift (as shown in Figs.|3|El and|SJ depends on 
the fact that states at large negative (positive) Eg cannot 
distort effectively due to having large numbers of satis- 
fied (unsatisfied) bonds that tend to have similar lengths. 
A system capable of uniform compression or expansion 
could take advantage of these states with extreme values 
of Eq in an entirely different manner than a system where 
the volume is constant. 
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